%{
AUTHOR: Felipe Arteaga

DATE:  2020--
-------------------------------------------------------------------------
PROJECT:
-------------------------------------------------------------------------
DESCRIPTION:

=========================================================================
%}

clc;clear;close all;fclose('all');
feature('DefaultCharacterSet', 'UTF-8');

pcName=char(java.lang.System.getProperty('user.name'));
if(strcmp(pcName,'felipe'))
    % PC Felipe
    myDir='/Users/felipe/Dropbox/';
    projectDir=[myDir,'git/warnings/'];
    projectDirData=[myDir,'projects/warnings/'];
end


addpath(genpath([myDir,'/myMatlabFunctions/']));

%%
dirPlots=[projectDir,'/paper/figuresCL/RDs/'];

% warning('Alternative dir plots')
% dirPlots=[projectDir,'/paper/figuresCL/RDs/Aux/'];

dirPlotsR=[projectDir,'/paper/figuresCL/RDs_multBandwidth/'];
dirTable=[projectDir,'/paper/tablesCL/'];

dirData=[projectDirData,'/data/chile/'];

heter=''; %'' , 'vulnerable', 'maturity ', 'sex'  '': none;
closePlotsWhilePlotting=true;

plotRDs=false;
savePlotRDs=false;
plotWithZoom=false;
bandwidthTable='local'; % 'full','local','localAuto';

     tableBandwidthComparison=false;
 

switch bandwidthTable
    case 'full'
        bwnote='Full bandwidth was used (-0.28,+0.68) with 2nd order (local) polynomial fit. ';
    case 'local'
        bwnote='Bandwidth was set to +-0.1 with default local polynomial fit (1st order). ';
    case 'localAuto'
        bwnote='Bandwidth was set to automatic (bwselect with default options) with default polynomial fit (1st order). ';
end

%% Load data:
anho=0; %0: pooled
if(anho==0)
    anhoStr='';
else
    anhoStr=sprintf('%i',anho);
end


notIn2020={}; % {'assignedToPref','acceptOffer','declineOffer','preferenciaAsign_1','placedInAnyPreferenceIniEnd','placedInAddedIniEnd','placedInAddedNoWaitlistIniEnd','placedInOldIniEnd','enrolledInAssigned','enrolledInAddedIniEnd','enrolledInAddedNoWaitlistIniEnd','enrolledInAddedNoRiskIniEnd'};

if(anho==0)
    
    d1=load([dirData,'2018/inputRD'],'dataRD');
    d2=load([dirData,'2019/inputRD'],'dataRD');
    d3=load([dirData,'2020/inputRD'],'dataRD');
    % Variables that are not available for 2020 yet:
    fillWithNan=notIn2020;
    %fillWithNan=
    for v=1:length(fillWithNan)
        d3.dataRD.(fillWithNan{v})=nan(height(d3.dataRD),1);
    end
    
    
    % Universe that is comparable between 2018-2020
    entryGrades=[-1 0 1 7 9];
    sirven1=ismember(d1.dataRD.grade,entryGrades)&not(d1.dataRD.cod_reg==13);
    sirven2=ismember(d2.dataRD.grade,entryGrades)&not(d2.dataRD.cod_reg==13);
    sirven3=ismember(d3.dataRD.grade,entryGrades)&not(d3.dataRD.cod_reg==13);
    
    d1.dataRD.comparable=sirven1;
    d2.dataRD.comparable=sirven2;
    d3.dataRD.comparable=sirven3;
    
    d1.dataRD.anho=2018*ones(height(d1.dataRD),1);
    d2.dataRD.anho=2019*ones(height(d2.dataRD),1);
    d3.dataRD.anho=2020*ones(height(d3.dataRD),1);
    
    % Keep vars in common
    
    incommon=intersect(intersect(d1.dataRD.Properties.VariableNames,d2.dataRD.Properties.VariableNames),d3.dataRD.Properties.VariableNames);
    incommon=incommon(not(strcmp(incommon,'mrun')));
    dataRD=[d1.dataRD(:,incommon);d2.dataRD(:,incommon);d3.dataRD(:,incommon)];
    
    
    
else
    load([dirData,anhoStr,'/inputRD'],'dataRD')
end

dataRD.rural=dataRD.newMarket<0;

% As asignment at .3 is control, just move (if any) in .3 to epsilon to the
% left:
dataRD.riskPopup(dataRD.riskPopup==.3)=.29999999;

% Avoid mass points at extremes:
dataRD.restrAll=dataRD.pobPopup==1&dataRD.riskPopup>.01&dataRD.riskPopup<.99;

% No need to be more specific, and put "Predicted placement risk of 1st
% attempt"
dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
dataRD.Properties.VariableDescriptions{'treatedPopup'}='Treated pop-up in first attempt';


if(false)
   %% % Popups of total population
   n_apps=[274990,483070,454415];
   
   fprintf('2018 In data: %.1f, in pop-up population: %.1f\n',height(d1.dataRD)/n_apps(1)*100,sum(d1.dataRD.pobPopup)/n_apps(1)*100);
   fprintf('2019 In data: %.1f, in pop-up population: %.1f\n',height(d2.dataRD)/n_apps(2)*100,sum(d2.dataRD.pobPopup)/n_apps(2)*100);
   fprintf('2020 In data: %.1f, in pop-up population: %.1f\n',height(d3.dataRD)/n_apps(3)*100,sum(d3.dataRD.pobPopup)/n_apps(3)*100);
   
   fprintf('\n ALL  In data: %.1f, in pop-up population: %.1f\n',height(dataRD)/sum(n_apps)*100,sum(dataRD.pobPopup)/sum(n_apps)*100);
   
end
   
   
  
% Add lista to enrolled rbd
for a=2018:2020
    dataLlenos=readtable(sprintf('%s/Mineduc/dataCompilada/Aux/_csv/llenos%i.csv',myDir,a));
    dataLlenos.anho=scalarForTable(a,dataLlenos);

    if(a>2018)
        dataLL=[dataLL;dataLlenos];
    else
        dataLL=[dataLlenos];
    end
end

convCodCursoCodNivel=readtable([myDir,'Mineduc/modelacion/riesgo/data/','2021','/otherData/relations/eqCodCursoACodNivelFake2019.csv'],'delimiter',',');
dataLL=outerjoin(dataLL,convCodCursoCodNivel,'keys',{'cod_curso'},'mergeKeys',true,'type','left','rightVariables',{'cod_nivel'});
dataLL.grade=dataLL.cod_nivel;
dataLL.rbdEnrolled=dataLL.rbd;


dataLLaux=stataCollapse({'rbdEnrolled','grade','anho'},dataLL,{'lleno','conLista'},{'mean'});
dataLLaux.oversubscribed=dataLLaux.conLista>.5;

assert(allunique(dataLLaux(:,{'rbdEnrolled','anho','grade'})))

dataRD=outerjoin(dataRD,dataLLaux,'keys',{'rbdEnrolled','anho','grade'},'mergeKeys',true,'type','left','rightVariables',{'oversubscribed'});
%dataRD.oversubscribed=scalarForTable(false,dataRD); 


%% RD

% Defines subpopulation in which I want to calculate RDs for outcomes
% defined soon.

if(anho==0) % This calculates estimates pooling all years
    
  
    
    switch heter
        case ''
            subpobs=cell(1,4);
            s=1;
            subpobs(s,:)={'All',    dataRD.restrAll,'pooled','addAnyIniEnd'};s=s+1;
            %subpobs(s,:)={'2018',   dataRD.restrAll&dataRD.anho==2018,'2018',''};s=s+1;
            %subpobs(s,:)={'2019',      dataRD.restrAll&dataRD.anho==2019,'2019',''};s=s+1;
            %subpobs(s,:)={'2020',    dataRD.restrAll&dataRD.anho==2020,'2020',''};s=s+1;
            posAll=1;
        
            
        otherwise
            error('aca')
    end
    
   
    
else
    subpobs=cell(1,3);
    s=1;
    subpobs(s,:)={'All',            dataRD.restrAll,'all',''};s=s+1;
    
end



% Definition of outcomes for RDs. 3rd row defines the beginning of a panel
% in the table
% 4th row

if(ismember('DRiskLExpostIniEnd',dataRD.Properties.VariableNames))
    error('Borrar esto');
else
    dataRD.DRiskExpostIniEnd=dataRD.riskExpostEnd-dataRD.riskExpostIni;
end

% New: enrolled conditional on placed (short name because of stata)
dataRD.enrolledInAssignedCOA=double(dataRD.enrolledInAssigned);
dataRD.enrolledInAssignedCOA(dataRD.assignedToPref==0)=nan;
dataRD.placedInAnyPreferenceNWIE=dataRD.placedInAnyPreferenceNoWaitlistIniEnd;
dataRD.placedInAnyPreferenceNoRiskIE=dataRD.placedInAnyPreferenceNoRiskIniEnd;



dataRD.valueAddedEnrolled(dataRD.grade>8)=nan;
dataRD.hasValueAddedEnrolled=double(not(isnan(dataRD.valueAddedEnrolled)));
dataRD.hasValueAddedEnrolled(dataRD.grade>8)=nan;

dataRD.valueAddedSep=dataRD.valueAddedEnrolled;
dataRD.valueAddedSep(not(dataRD.prioritario))=nan;
dataRD.valueAddedPlaced=dataRD.valueAddedEnrolled;
dataRD.valueAddedPlaced(not(dataRD.placedInAnyPreferenceIniEnd))=nan;
dataRD.valueAddedNotPlaced=dataRD.valueAddedEnrolled;
dataRD.valueAddedNotPlaced((dataRD.placedInAnyPreferenceIniEnd))=nan;



dataRD.undersubscribed=not(dataRD.oversubscribed);


dataRD.overAndEnrolledNotHS=double(dataRD.oversubscribed);
dataRD.overAndEnrolledNotHS(isnan(dataRD.rbdAnhoSig))=nan;
dataRD.overAndEnrolledNotHS(isnan(dataRD.grade>8))=nan;

dataRD.overAndEnrolledNotHSwithVA=dataRD.overAndEnrolledNotHS;
dataRD.overAndEnrolledNotHSwithVA(isnan(dataRD.valueAddedEnrolled))=nan;

dataRD.residualVA=scalarForTable(nan,dataRD);
dataRD.residualVA(dataRD.undersubscribed)=dataRD.valueAddedEnrolled(dataRD.undersubscribed)-mean(dataRD.valueAddedEnrolled(dataRD.undersubscribed),'omitnan');
dataRD.residualVA(dataRD.oversubscribed)=dataRD.valueAddedEnrolled(dataRD.oversubscribed)-mean(dataRD.valueAddedEnrolled(dataRD.oversubscribed),'omitnan');

dataRD.expectedVA=scalarForTable(nan,dataRD);
dataRD.expectedVA(dataRD.undersubscribed)=mean(dataRD.valueAddedEnrolled(dataRD.undersubscribed),'omitnan');
dataRD.expectedVA(dataRD.oversubscribed)=mean(dataRD.valueAddedEnrolled(dataRD.oversubscribed),'omitnan');
dataRD.expectedVA(isnan(dataRD.valueAddedEnrolled))=nan;



dataRD.residualVA2=scalarForTable(nan,dataRD);
dataRD.residualVA2(dataRD.undersubscribed)=dataRD.valueAddedEnrolled(dataRD.undersubscribed)-mean(dataRD.valueAddedEnrolled(dataRD.undersubscribed&dataRD.restrAll),'omitnan');
dataRD.residualVA2(dataRD.oversubscribed)=dataRD.valueAddedEnrolled(dataRD.oversubscribed)-mean(dataRD.valueAddedEnrolled(dataRD.oversubscribed&dataRD.restrAll),'omitnan');

dataRD.residualVA_under=scalarForTable(nan,dataRD);
dataRD.residualVA_under(dataRD.undersubscribed)=dataRD.valueAddedEnrolled(dataRD.undersubscribed)-mean(dataRD.valueAddedEnrolled(dataRD.undersubscribed),'omitnan');

dataRD.residualVA_over=scalarForTable(nan,dataRD);
dataRD.residualVA_over(dataRD.oversubscribed)=dataRD.valueAddedEnrolled(dataRD.oversubscribed)-mean(dataRD.valueAddedEnrolled(dataRD.oversubscribed),'omitnan');

%Y_j=E[Y_j|type]+1[type=over]*(Y_j-E[Y_j|type=over)+1[type=under]*(Y_j-E[Y_j|type=under);

% 5th column defines if IV is calculated or not.
outcomes={'addAnyIniEnd','Add any','aa','',0;...
    'valueAddedEnrolled','Value added (VA)','v1','',1;...
    'overAndEnrolledNotHS','Enrrolled in oversubscribed ($type=over$)','v1_5','',1;...
    'overAndEnrolledNotHSwithVA','Enrrolled in oversubscribed ($type=over$)|not missing VA','v1_6','',1;...
    'expectedVA','$E[VA|type]$','v0','',1;...
    'residualVA','$VA-E[VA|type]$','v2','',1;...
    %'residualVA2','Residual2','v3','',0;...
       'residualVA_under','$(VA-E[VA|type])|type=under$','v4','',1;...
          'residualVA_over','$(VA-E[VA|type])|type=over$','v5','',1;...
    };


% No yeaar and IV.




% outcomes={   % 'riskLastDayIni','True Risk initial attempt','ri','';...
%     %'riskLastDayEnd','True Risk final attempt','rf','';...
%     %'DRiskLastDayIniEnd','$\Delta$ True Risk','dr','';...
%     'riskExpostIni','True Risk initial attempt','riEP','';...
%     'riskExpostEnd','True Risk final attempt','rfEF','';...
%     'DRiskExpostIniEnd','$\Delta$ True Risk ExP','drEP','';...
%     'placedInAnyPreferenceIniEnd','Placed to preference','ap','C. Choice outcome';...
%     'placedInAddedIniEnd','Placed to added preference','apad','';...
%     'placedInOldIniEnd','Placed to old preference','apold','';...
%     };


if(anho==2020)
    outcomes=outcomes(not(ismember(outcomes(:,1),notIn2020)),:);
end



dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
for o=1:size(outcomes,1)
    
    dataRD.Properties.VariableDescriptions{outcomes{o,1}}=outcomes{o,2};
end


assert(allunique(outcomes(:,1)))
assert(allunique(outcomes(:,2)))
assert(allunique(outcomes(:,3)))


data=dataRD(:,[{'riskPopup'},outcomes(:,1)']);
inAnyPob=false(height(data),1);


% Generate stata commands
commands=cell(size(subpobs,1)*size(outcomes,1),2);
counter=1;

for p=1:size(subpobs,1)
    
    data.(sprintf('subpop%i',p))=subpobs{p,2};
    inAnyPob=inAnyPob|subpobs{p,2};
    
    for o=1:size(outcomes,1)
        
        % Avoid computing outcomes that are not available:
        noCalcular=(strcmp(subpobs{p,3},'2020')|strcmp(subpobs{p,3},'2020comp'))&ismember(outcomes{o,1},notIn2020);
        
        if(not(noCalcular))
       
            % Full bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'full')||plotRDs)
                commands(counter,:)={sprintf('%s_%i_full',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, c(.3) p(2) h(.28 .68)',outcomes{o,1},p)};counter=counter+1;
            end
            % Full loccal fixed bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'local'))
                commands(counter,:)={sprintf('%s_%i_local',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, c(.3)  h(.1 .1)',outcomes{o,1},p)};counter=counter+1;
            end
            % Full local auto bandwidht
            if(tableBandwidthComparison||strcmp(bandwidthTable,'localAuto'))
                commands(counter,:)={sprintf('%s_%i_localAuto',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, c(.3) ',outcomes{o,1},p)};counter=counter+1;
            end
            
            % If with IV:
            if(outcomes{o,5}==1&&~isempty(subpobs{p,4}))
                
                endogenousVar=subpobs{p,4};
                
                % Full bandwidth
                if(strcmp(bandwidthTable,'full'))
                    commands(counter,:)={sprintf('%s_%i_full_iv',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, fuzzy(%s) c(.3) p(2) h(.28 .68)',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
                % Full loccal fixed bandwidth
                if(strcmp(bandwidthTable,'local'))
                    commands(counter,:)={sprintf('%s_%i_local_iv',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, fuzzy(%s) c(.3)  h(.1 .1)',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
                % Full local auto bandwidht
                if(strcmp(bandwidthTable,'localAuto'))
                    commands(counter,:)={sprintf('%s_%i_localAuto_iv',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, fuzzy(%s) c(.3) ',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
            end
        end
    end    
end
commands=commands(not(cellfun(@isempty,commands(:,1),'UniformOutput',true)),:);


if(false)
    % Run specific rd
    %%
dataRD.auxRest=dataRD.restrAll;
a=stataCommand('rdrobust valueAddedEnrolled riskPopup if auxRest==1, c(.3) p(1) h(.1 .1)',dataRD);
a.log
%stataExpress('rdrobust  DRiskExpostIniEnd riskPopup if auxRest==1, fuzzy(addAnyIniEnd) c(.3) p(1) h(.1 .1)',data)
end

% Run Stata:
data=data(inAnyPob,:);
res=stataCommand(commands,data);


%% Make plots and table


cantPobs=(size(subpobs,1));
cantRegs=(size(outcomes,1));

matrix_nl=nan(cantRegs,cantPobs);
matrix_nr=nan(cantRegs,cantPobs);
matrix_b=nan(cantRegs,cantPobs);
matrix_se=nan(cantRegs,cantPobs);

matrix_b_iv=nan(cantRegs,cantPobs);
matrix_se_iv=nan(cantRegs,cantPobs);
matrix_nl_iv=nan(cantRegs,cantPobs);
matrix_nr_iv=nan(cantRegs,cantPobs);

matrix_b_l=nan(cantRegs,cantPobs);
matrix_se_l=nan(cantRegs,cantPobs);

matrix_b_r=nan(cantRegs,cantPobs);
matrix_se_r=nan(cantRegs,cantPobs);

matrix_biasCorr_ci=nan(cantRegs,cantPobs,2);



regNames=fieldnames(res);

for p=1:cantPobs
    for r=1:cantRegs
        
        regName_i=sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable);
        if(ismember(regName_i,regNames))
            res_i=res.(regName_i);
            
            % Beta
            matrix_b(r,p)=res_i.tau_cl;
            % Sd Beta
            matrix_se(r,p)=res_i.se_tau_cl;
            
            % Beta left
            matrix_b_l(r,p)=res_i.beta_p_l(1);
            % Sd Beta left
            matrix_se_l(r,p)=sqrt(res_i.V_cl_l(1));
            
            % Beta right
            matrix_b_r(r,p)=res_i.beta_p_r(1);
            % Sd Beta right
            matrix_se_r(r,p)=sqrt(res_i.V_cl_r(1));
            
            % Bias corrected Confidence interval (alpha=5%)
            matrix_biasCorr_ci(r,p,1)=res_i.tau_bc-norminv(.975)*res_i.se_tau_rb;
            matrix_biasCorr_ci(r,p,2)=res_i.tau_bc+norminv(.975)*res_i.se_tau_rb;
            
            % N
            matrix_nl(r,p)=res_i.N_h_l;
            matrix_nr(r,p)=res_i.N_h_r;
            
            % If with IV:
            if(outcomes{r,5}==1&&~isempty(subpobs{p,4}))
                regName_i_iv=sprintf('%s_%i_%s_iv',outcomes{r,3},p,bandwidthTable);
                
                res_i_iv=res.(regName_i_iv);
                % Beta
                matrix_b_iv(r,p)=res_i_iv.tau_cl;
                % Sd Beta
                matrix_se_iv(r,p)=res_i_iv.se_tau_cl;
                
                % N
                matrix_nl_iv(r,p)=res_i_iv.N_h_l;
                matrix_nr_iv(r,p)=res_i_iv.N_h_r;
                
            end
            
            if(plotRDs)
                
                % Load full bandwidth result for ploting
                warning('Add small bandwidth estimate for full bandwidth plot')
                
                res_i=res.(sprintf('%s_%i_%s',outcomes{r,3},p,'full'));
                res_i_table=res.(sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable));
                
                plotsub=subpobs{p,2};
                
                figure
                plotRD(res_i,dataRD,plotsub,'otherPointEstimate',res_i_table);
                if(savePlotRDs)
                    easyExport([dirPlots,sprintf(sprintf('%s_P_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                
                if(plotWithZoom)
                    % Plot with zoom
                    figure
                    plotRD(res_i,dataRD,plotsub,'newxlim',[.1,.5],'nb',100,'otherPointEstimate',res_i_table);
                    easyExport([dirPlots,sprintf(sprintf('%s_P_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                if(closePlotsWhilePlotting)
                    close all
                end
            else
                fprintf('Ojo: reg %s no existe para subpob %s\n',regName_i,subpobs{p,1});
            end
        end
    end
    
    %vector_Nl(1,p)=sum(dataRD.riskPopup<.3&subpobs{p,2});
    %vector_Nr(1,p)=sum(dataRD.riskPopup>.3&subpobs{p,2});
    
end

for p=1:cantPobs
    %vector_Nl(1,p)=sum(dataRD.riskPopup<.3&subpobs{p,2});
    %vector_Nr(1,p)=sum(dataRD.riskPopup>.3&subpobs{p,2});
end




%% RD main table
matTable=nan(size(matrix_b).*[1 2]+[2 0]);
matTable_sd=nan(size(matrix_b).*[1 2]+[2 0]);

matTable(:,1:2:end)=[matrix_b;...
    matrix_nl(2,:)...
    ;matrix_nr(2,:)];

if(any(not(isnan(matrix_b_iv)),'all'))
matTable(:,2:2:end)=[matrix_b_iv;...
    matrix_nl_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)...
    ;matrix_nr_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)];
end

matTable_sd(:,1:2:end)=[matrix_se;nan(2,size(matrix_se,2))];
matTable_sd(:,2:2:end)=[matrix_se_iv;nan(2,size(matrix_se_iv,2))];

withInfo=not(all(isnan(matTable),1));


header=repmat({''},2,size(matTable,2));
header(1,1:2:end)=subpobs(:,1)';
header(1,2:2:end)=subpobs(:,1)';
header(2,2:2:end)={'IV'};

matTable=matTable(:,withInfo);
matTable_sd=matTable_sd(:,withInfo);
header=header(:,withInfo);


opt=struct;
opt.header=header;
opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=[outcomes(:,2);{'NL';'NR'}];
%opt.stars=getStars(matTable,matTable_sd);
opt.addColumnNumber=true;
if(anho==0)
    %opt.columnafantasma=4;

else
    opt.title=sprintf('%s Sample - Pop-Up effect',anhoStr);
end
opt.filaFantasma=size(matTable,1)-2;
paneles=find(not(ismissing(outcomes(:,end-1))));
opt.panel=[num2cell(paneles-1),outcomes(paneles,end-1)];
%opt.alignmentFirstCol={'L{3cm}'};
opt.adjust=true;

%opt.file=sprintf('%s%s_tablePopup%s_play',dirTable,anhoStr,heter);

%        opt.columnafantasma=4;
        opt.title='RD estimates of platform pop-up effects';
        %opt.label='tab:RDpopup';
        opt.note=['Local linear RD estimates of pop-up effects from warning pop-up on application platform. Computed using triangular kernel with bandwidth 0.1. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. ',...
            'We report estimates in the pooled sample across years 2018-2020. IV (column 2) shows the instrumental variable specifications, where the endogenous regressor is the add any school indicator. ',...
        ];
tabla=cell2latex(mat2cellstr(matTable,'precisionDecimal','%.3f','revisarFilas',true),'opts',opt);
compileLatex(tabla)


